NASA REPORT NO. CR 174926 


ylRVIN/C/H-SFVIN 


A HYBRID NUMERICAL TECHNIQUE FOR 
PREDICTING THE A EROD YNA MIC A ND 
ACOUSTIC FIELDS OF ADVANCED TURBOPROPS 

Contract No. NAS3-23699 


April 1985 


FINAL REPORT 

CALSPAN REPORT NO. 7157-A-1 


Prepared by: 

G.F. Homicz and J.R. Moselle 
PHYSICAL SCIENCES DEPARTMENT 


Prepared for: 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
LEWIS RESEARCH CENTER 
CLEVELAND, OHIO 44135 


CALSPAN CORPORATION 


P.O. BOX 400 BUFFALO. NEW YORK 1 4225 


FOREWORD 


This is the final report documenting the results of a theoretical 
program to predict the aerodynamic and acoustic performance of advanced 
turboprops, sponsored by the NASA Lewis Research Center. The program extended 
from April 1983 to March 1985, with Dr. Kenneth Baumeister serving as Technical 
Monitor. 


Dr. Paul V. Marrone, Head of the Physical Sciences Department, has 
overall responsibility for management and review of all technical programs 
within the department. It was originally envisioned that Dr. William J. Rae 
would serve as Principal Investigator, with assistance from Drs. Gregory F. 
Homicz, Joseph P. Nenni, and John A. Lordi. Shortly after the program began, 
the departure of Drs. Rae, Nenni and Lordi led to Dr. Homicz being appointed 
Principal Investigator. To aid him in the technical effort, Calspan obtained 
the services of Prof. A. Seybert of the University of Kentucky as a 
subcontractor. Prof. Seybert had responsibility for the outer linearized 
acoustic analysis and the preparation of Section 4 and portions of Section 5 of 
this report. We wish to thank Dr. L. Bober, M. Celestina, and H. Huynh of NASA 
Lewis for their generous help in familiarizing us with the NASPROP-E code and 
their system. 
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ABSTRACT 


A hybrid numerical procedure is presented for the prediction of the 
aerodynamic and acoustic performance of advanced turboprops. Because the 
relative tip speed at design is typically supersonic, one must anticipate the 
presence of shock waves as well as transonic three-dimensional effects in the 
immediate vicinity of the blades. The strong rotational character of such flows 
warrants the use of the inviscid nonlinear Euler equations in predicting 
aerodynamic performance. In the outer flow regime away from the propeller such 
effects have decayed significantly; here the primary concern is the sound being 
propagated to the farfield. Hence in this region a linearized acoustic analysis 
is justified. 

The present investigation proposes a hybrid scheme which in principle 
leads to a consistent simultaneous prediction of both fields. In the inner flow 
a finite difference method, the Approximate-Factorization Alternating-Direction- 
Implicit (ADI) scheme, is used to solve the nonlinear Euler equations. In the 
outer flow the linearized acoustic equations are solved via a Boundary-Integral 
Equation (BIE) method. The two solutions are iteratively matched across a 
fictitious interface in the flow so as to maintain continuity. At convergence 
the resulting aerodynamic load predictions will automatically satisfy the 
appropriate "free-field" boundary conditions at the edge of the finite 
difference grid, while the acoustic predictions will reflect the "back -reaction" 
of the radiated field on the magnitude of the loading source terms, as well as 
refractive effects in the inner flow. 

The equations and logic needed to match the two solutions are 
developed and the computer program implementing the procedure is described. 
Unfortunately, no converged solutions were obtained, due to the unexpectedly 
large running times. The reasons for this are discussed and several means to 
alleviate the situation are suggested. 
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Section 1 
INTRODUCTION 


The past two or three decades have witnessed a progression in the 
commercial aviation fleet from propeller to turbojet, and most recently, to 
turbofan-powered aircraft. The transition has been motivated primarily by the 
desire for increased speed, reduced cabin vibration, and relatively simpler 
design. Early turbojet engines were rather noisy, due to the very high exit 
velocities of the turbulent exhaust jet. This eventually led to the turbofan 
engine, whose exhaust velocity (and resultant noise) is much lower. To 
compensate for the reduced jet thrust, a significant fraction of the inlet air 
is made to bypass the core engine and flows instead through a ducted fan, which 
generates a significant portion of the thrust. 

In recent years the need for increased fuel efficiency has prompted 
reconsideration of the propeller as a viable alternative. Operation at a 
typical cruise condition (M = 0.8 at 30,000 ft.) usually means blade relative 
tip speeds which are in the transonic/low supersonic regime. The performance 
degradation and "buzz-saw" noise produced by conventional turboprops at these 
conditions preclude their use. Recent advances in our understanding of such 
flows, and improvements in manufacturing techniques, have driven propeller 
design in the direction of using many smaller diameter, thinner blades, with a 
highly swept planform. The similarity to an unducted fan is obvious, hence the 
name "propfan". 

While wind-tunnel as well as flight tests on scale models have 
confirmed the propulsive efficiency of modern propfan design, the noise they 
create is still a major concern. Most noise prediction methodologies model the 
blades as distributions of source and dipole acoustic singularities. The 
magnitude of these singularities is determined by the blades' thickness and 
loading distributions, respectively. In principle the geometric thickness 
distribution is a given, but the loads themselves require another whole 
prediction methodology, which for advanced propfans is generally a numerical 
algorithm of some sort. Usually the numerical solution is extended only a 
finite distance from the propeller, where uniform free-stream conditions are 
assumed to exist. This assumption is clearly inconsistent with acoustic 
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radiation to the farfield. Such an approach thus completely neglects any 
coupling between the aerodynamic and acoustic field, i.e. the "back -react ion" of 
the radiative field on the source mechanism. 

Comparisons of such noise predictions against available wind-tunnel 
and flight test acoustic data, e.g. Refs. 1 and 2, show uneven agreement. The 
Sound Pressure Level (SPL) for both theory and experiment exhibits a nearly 
linear increase with increasing tip speed in the subsonic regime. However, as 
the relative tip Mach number passes through unity, the data tend to level off at 
a plateau, while the predicted SPL continues to rise. It is the discrepancy 
between theory and experiment in the transonic regime which this program is 
intended to address. Specifically, we hope to determine whether a direct 
coupling of the aerodynamic and acoustic fields in the theory would improve the 
agreement . 


For propeller flows which are subsonic at the hub and supersonic at 
the tip, strong transonic and nonlinear effects must be anticipated, including 
three-dimensionality, swirl, and even shock waves. These phenomena dictate the 
solution of the three-dimensional, nonlinear, inviscid Euler equations if 
accurate aerodynamic predictions are required. For this purpose we have chosen 
as our starting point the NASPROP-E computer code written by Chaussee and Kutler 
(Ref. 3)» This choice was dictated by the code's ready availability on the NASA 
Lewis CRAY system, lack of proprietary restraints on its use or modification, and 
its demonstrated success in treating advanced propfan configurations (Ref. 4). 

NASPROP-E solves the Euler equations using a finite-difference 
algorithm on a boundary-conforming grid. In principle the grid and solution 
could be extended to the acoustic farfield, but in practice this would soon 
stretch storage requirements and CPU times to the breaking point. Out of 
necessity the outer boundary of the grid is placed a finite distance, say 
several blade diameters, away from the propeller where the boundary conditions 
are only approximated by undisturbed free-stream values. 

On the other hand, one expects intuitively that the nonlinear effects 
necessitating the use of the full Euler equations are likely to be dominant 
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only in the immediate vicinity of the blades. In the farfield, where the 
propagation of acoustic waves is the prime concern, an analysis based on the 
linearized flow equations should suffice. The latter lend themselves to 
semi -analytical methods which have a less voracious appetite for computer 
resources. This suggests a hybrid methodology in which NASPROP-E is used in the 
inner flow and matched across some suitably chosen interface to a linearized 
model of the outer flow. 

The scheme chosen here for the outer flow is the Boundary-Integral 
Equation (BIE) method (Ref. 5). As the name suggests, it is based on converting 
the partial differential equations to an equivalent integral equation on the 
interface. The number of spatial dimensions is thus reduced from three to two. 
As envisioned here, the NASPROP-E solution in the inner flow is used to compute 
the pressure gradient 9-p/dn (n being normal to the interface). This is used 
as the inner boundary condition for the outer flow solution. The latter will 
uniquely determine a new p distribution on the interface, fully consistent with 
the appropriate radiation condition ab infinity. This new p distribution is 
then used as the outer boundary condition on a new calculation of the inner 
flow. This procedure is repeated until convergence is reached, as evidenced by 
continuity of pressure and its normal derivative across the interface. After 
convergence, the SPL at any point in the outer flow is available as a simple 
quadrature of p and df>/dn (which are now known) over the interface. 

The discussion in the following sections proceeds logically from the 
propeller outwards. Section 2 gives an overview of the equations, variables, 
and algorithm used in NASPROP-E. Section 3 discusses how this solution is 
matched across the interface to the BIE solution, which in turn is discussed in 
Section 4. The overall structure of the combined code, termed CALPROP, is 
presented in Section 5 along with a description of the additional inputs and 

outputs. Large computer run times have thwarted efforts to get a converged 
solution, the reasons for which are discussed in Section 6. Section 7 closes 
with the major conclusions from this study and suggestions for future research. 
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Section 2 

INNER NONLINEAR AERODYNAMIC MODEL 

As noted in the Introduction, advanced propfan configurations can be 
expected to exhibit strong rotational flow effects, as well as shock waves. For 
this reason a solution of the full nonlinear, inviscid, three-dimensional Euler 
equations is needed in the immediate vicinity of the blades. As our starting 
point, we chose the NASPROP-E code developed by Chaussee and Sutler (Ref. 3)* 
This source code was readily available and working on the NASA Lewis CRAY 
system, and has been successfully used to analyze configurations of interest 
here (Ref. 4). 

Actually, the overall package is a collection of three stand-alone 
programs, each designed to do a different task in the analysis sequence. The 
first is a mesh generation program which takes the propfan/nacelle geometry as 
input and produces a suitable boundary-conforming grid. This grid is used as 
input to the flowfield analysis program, which solves the equations and stores 
the full three-dimensional solution at all grid points onto a disk file. The 
last step is the data reduction program which reads the solution on file and 
computes force and moment coefficients, generates plots of selected quantities, 
etc. The flowfield analysis program is by far the most complex and time 
consuming. It is this portion that we are primarily concerned with here, with 
the exception of one small change noted below in connection with the mesh 
generation code. 

As NASPROP-E has been documented elsewhere, we give here only an 
overview sufficient to introduce certain variables and concepts. We also 
discusss what modifications were found necessary to allow interaction with the 
outer* solution field. 

Grid Transformation 


It is assumed that the inlet flow, at subsonic Mach number, M, is 
uniform. Hence to an observer in a blade-fixed reference frame (2,?~, <fi ) the 
flow appears steady, and it is in such a frame that the equations are solved. 
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Further, since all of the B blades are assumed identical, it is only necessary 
to solve for the flow through one blade passage , say O ^ (f> . To 

facilitate implementation of the blade and nacelle boundary conditions, the flow 
is mapped into computational coordinates, 

4 - 4 

K a h (2,r,<p, t) (2-1) 

C *) 

't • t 

The transformation is boundary-conforming, i.e., the flow boundaries in the 
computational domain all lie along constant values of one of the coordinates. 

In this domain the flow occupies a rectangular prism, within which a uniform 
Cartesian grid is employed. Such a transformation makes it easy to adapt any 
number of finite-difference algorithms to the solution. It also allows one* to 
"bunch up" grid points in regions of strong gradients without compromising the 
accuracy of the difference scheme. 

• A cross-section in the ( Z , ) plane of the difference grid used 

here is shown in Fig. 1. This displays the non-uniform curvilinear grid in 
physical space corresponding to a uniform Cartesian grid in <£ and >^. The 
corresponding integer indices are J and K, respectively. J = 1 represents the 
horizontal stagnation streamline coming in from the left. Higher values of J 
are the lines rotating to the right, with J : JA denoting the specific <4 = 
constant line intersecting the upper left corner of the grid. The vertical 
outflow boundary is denoted by J = JMAX. Note the bunching near the stagnation 
point and the blade leading and trailing edges . The intersecting grid lines are 
for >i and K constant; K = 1 is the nacelle surface, while K = KMAX along the 
vertical inflow boundary and the horizontal cylindrical sidewall. In the 
azimuthal direction, L is the index for coordinate C, (not shown in Fig. 1, see 
Ref. 4); L = 1 lies along the blades' suction surface, and L = LMAX along the 
blades' pressure surface. 

Transformed Equations 

In the computational space the Euler equations can still be written in 
weak -conservation law form as follows: 


KMAX 
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Figure 1. FINITE DIFFERENCE GRID IN z, r PLANE 




where subscripts are used to denote differentiation. Q is a column vector 
consisting of the unknowns , 

Q - [j> t J>tL, eJ T ( 2 - 3 ) 


where is the fluid density, (u, v, w) the velocity components in the 
directions, respectively, and e is the total energy per unit volume defined by 


e 


P 

Cr-0 


+ 


JO (u. z + zr z + 1<S Z ) 
_ 


(2-4) 


where p is the pressure. The quantity J is the Jacobian of the transformation 
defined in Eq. (2-1), i.e. 




and Y is the ratio of specific heats for an ideal gas . 


(2-5) 


Before going further, it should be noted that all variables in the 
code have been nondimensionalized. Lengths are normalized by the propeller 
diameter, D. Pressure and density are normalized by their freestream values, p^ 
and . Velocities are normalized by 0. eo j~\fY' > where Q.^ is the free-stream 
sound speed, • Time is normalized by DY Y'f . Finally, from 

Eq. (2-4) the energy per unit volume is also normalized by p^ . 

The quantities E, F and G are the column flux vectors in the <£ , y[ and £ 
directions defined in terms of the unknowns as , 

E = J 'r/>V, puU+p £ z , p , j2url/+ Jbtf //-, (e+p)CS'p£ t ] T 

F = S'* [p\/ t />uV+p^, fv-\Z+pn r ,purV+f3fy/r,(e+-p)V-p>i t ] T 

G = [p b/ } pul , fv-lV+pS r) ft*slV+f<;t/r, (e + f>) £J r 

( 2 - 6 ) 

Here (U, V, W) are the contravariant velocity components in the ( £ , , <£ ) 

directions, respectively: 
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(2-7) 


U - 4+ * U< ?B + 

V = y[ t + “fa + v >i r + 

W - Cj. * u£ m + v $r * 

The undifferentiated source term H is given by 

H » ('J' fr )~[fV, jOUV-jjO (v x -Us z ) } z/>vzo-j (e+f>jzrj ( 2 - 8 ) 

The metrics of the transformation appearing in Eqs. (2-6) and (2-7) can also be 
evaluated in the computational space using the chain rule, e.g. 

J f^0c ' < P^) 

Z r ~ ~ z \4s) C2-9) 

i* - j ( z>t+k 

and similarly for the >7 and £ metrics. 

Numerical Algorithm 

The vector Eq. (2-2) is a shorthand for 5 scalar equations 
representing, respectively, conservation of mass, conservation of momentum (3 
components), and conservation of energy. Note that the time derivative has been 
retained even though we seek the steady-state solution. This renders the 
equations hyperbolic and thus amenable to a time-marching algorithm. That is, 
the solution is integrated from some (arbitrary) initial condition through a 
sequence of finite time steps until all transients have decayed. This 

asymptotic field is then the desired steady-state solution. 

The use of an explicit finite-difference algorithm, though 
straightforward, is hampered by the very small time steps required to maintain 
numerical stability when small grid spacings are present. This is particularly 
nettlesome in applications such as the present, where the transients are not of 
interest, and we wish to converge to the steady solution as quickly as possible. 
For this reason NASPROP-E employs an implicit marching algorithm. Though more 
operations are required at each time step, implicit schemes allow greatly 
increased step sizes without sacrificing stability. The net result is a 
significant savings in CPU time. 
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The particular scheme used in NASPROP-E is the Approximate- 
Factorization Alternating-Direction-Implicit (ADI) scheme in so-called "delta" 
form (Refs. 6, 7). Basically, after applying differencing formulae to Eq. 

(2-2), the resulting three-dimensional spatial difference operator is factored 
into the product of three one -dimensional operators. Each of these operators 
represents a block -tridiagonal matrix. The latter can be inverted one at a time 
sequentially, a significant savings in both storage and time requirements 
compared with inverting the much larger matrix for the full three-dimensional 
operator. A method for speeding the calculation even more through further 
diagonalization of the block matrix structure (Refs. 8, 9) is also used. For 
more details on the Approximate -Factorization ADI algorithm and its 
implementation in NASPROP-E, see Refs. 3» 4, 6-9. No changes to this part of 
the code were made during the present investigation. 

Initial and Boundary Conditions 

The usual initial condition when starting a case from scratch is to 
set all the flow variables in the vector Q equal to their free-stream values. 
Optionally, one can also start the calculation using the flowfield solution 
stored as the result of a previous run. No changes have been made to the code 
in this regard (Section 4.4 of Ref. 3)* 

The boundary conditions used are those appropriate to an inviscid 
calculation. Since the blades and nacelle are impermeable the surface tangency 
condition is invoked, i.e., the velocity normal to the surface must be zero. 

In the computational space the blades map to a surface of constant C > so that 
the corresponding contravariant velocity, W, must vanish. Similarly, the 
nacelle maps to a surface of constant , along which the contravariant velocity 
V must be zero. Details of how these are enforced may be found in Ref. 3; no 
changes were made for the present application. 

As originally written, the outer boundary of the NASPROP-E grid was 
not as depicted here in Fig. 1. The upstream portion, rather than terminating 
in a flat vertical surface, consisted of a hemispherical cap of the same radius 
as the downstream cylindrical portion. This will be referred to as the 
hemisphere-cylindrical grid. On the hemispherical cap and cylindrical sidewall, 
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K = KMAX , free-stream values for all the flow quantities were imposed for all 
time. On the vertical downstream outflow boundary (J = JMAX) only the pressure 
was assumed to have returned to a constant value of p M . Significant variations 
in the velocity field must be expected, however, because of the added momentum 
and swirl imparted to the flow. Appropriate values for the velocity components 
were computed using the Method of Characteristics (MOC); see Ref. 3 for details. 
The density was obtained from the isentropic relation, and the total energy per 
•unit volume then follows from Eq. (2-4). 

The above scheme, assuming as it does a constant value for the 
pressure over the whole boundary of the grid, is clearly inconsistent with 
acoustic radiation to the farfield. It is this portion of NASPROP-E which had 
to be modified for the present application. First, the shape of the outer 
boundary in physical space was changed from the hemisphere-cylinder combination 
described in the last paragraph, to the fully cylindrical shape shown in Fig. 1. 
Although the matching to the outer flow solution can in principle be 
accomplished for an arbitrarily shaped interface, for reasons which will become 
clearer in Section 3 it is much easier if each portion of the interface is 
parallel to one of the coordinate directions. We are indebted to Drs. L. Bober 
and H. Huynh of NASA Lewis who made the necessary mesh alterations for us and 
stored the new grid on their CRAY system where we could directly access it. It 
should be noted that this change involves only the mesh generation program, and 
is essentially transparent to the three-dimensional flowfield program. 

Changes have also been made within the flowfield program to allow for 
a variable pressure distribution on the interface. As explained in Sections 
3 and 4, NASPROP-E is used to compute d-pjdri on the interface, which is used by 
the BIE package to compute a new p distribution there. This new information 
must somehow then be used as an updated outer boundary condition for the next 
iteration through NASPROP-E; this is complicated by the fact that p is not one 
of the primary variables stored in the solution vector Q, cf. Eq. (2-3). 
Accordingly, at the upstream vertical face and the cylindrical sidewall in Fig. 

1 the updated pressure from BIE is used in Eq. (2-4) to compute and store new 
total energy values on the interface, assuming the velocity components retain 
their free-stream values. On the downstream vertical face, where, u, v and w 
vary significantly, the same MOC method is still used to set new values for the 
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velocities, density and energy. However, now the p values from BIE are used as 
input rather than p^ . 

This completes the discussion of the NASPROP-E portion of the code. 
The next section discusses in detail the theory behind its interaction with the 
BIE solution of the outer flow. 
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Section 3 

MATCHING ACROSS THE INTERFACE 

The preceding section gave an overview of how the NASPROP-E code 
solves the nonlinear Euler equations in the inner flow region, and what changes 
were necessary to allow it to accept arbitrary outer pressure boundary conditions 
from the BIE package. This section discusses in more detail the theory behind 
how the two codes were intended to interact. There were three major 
incompatibilities between the codes which had to be addressed: NASPROP-E 

originally used a hemisphere/cylinder outer boundary, while BIE was designed 
with a fully cylindrical inner boundary in mind; NASPROP-E works in blade-fixed 
coordinates whereas the coordinate system in BIE does not rotate; NASPROP-E 
allows for a uniform subsonic axial flow, while the classical acoustic analysis 
in BIE does not. The first of these was easily resolved by switching to a 
cylindrical outer boundary in NASPROP-E, as described in Section 2. As will be 
seen shortly, this also facilitates transfer of the boundary conditions between 
the two. Resolution of the remaining incompatibilities is described below. 

Evaluation of Normal Pressure Gradient at Interface 


Figure 2 is a sketch of the same cylindrical interface shown in Fig. 

1 , but here with emphasis on the outer flow region . S A refers to the upstream 
vertical face at axial location z u , Sg to the cylindrical sidewall of radius Rg, 
and S c to the downstream vertical face at a^; the total ensemble is denoted by 
S. Keep in mind that S is a fictitious surface sufficiently far removed from 
the propeller, several blade diameters say, so that nonlinear effects such as 
shock waves are assumed to be negligible. 

We have the choice of using NASPROP-E to specify d*p/dn on S and then 
using BIE to get p, or using NASPROP-E to specify p on S and using BIE to 
determine d*pfdn . Due to the way in which the outer boundary conditions are 
set in NASPROP-E (cf. Section 2), the first alternative is clearly preferable. 

We thus need to evaluate the outward normal pressure gradient in the blade -fixed 
( Z , ir , (p ) frame in terms of the solution computed in the boundary-conforming 
( d, , >1 ,C ) coordinates. Using the chain rule for differentiation on surface 
S A we have 

* 7* A * *•?*) (3-la) 
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ARE MATCHED (not to scale) 



while on S 3 , 
and on S c , 

df)/d/7 - = 

The metrics ,£ r etc. are already stored (cf. Eq. (2-9)). Evaluation of the 
pressure gradients on the right is a two-step process. First the pressures at 
grid points on and near S are evaluated for the stored solution vector Q using 
Eq. (2-4). Then standard second-order accurate difference formulae are applied 
to get p^ , and p^ . 

Equations (3-1) define the normal pressure gradient only within the 
pie-shaped sector formed by one blade passage, say O £ -y~$Q s and 
Because all the blades are identical, this pattern will repeat itself B times 
around the axis in both hub-fixed and blade-fixed coordinates. It is natural 
then to represent the variation as a Fourier series: 

*)/*» “ (3 ' 2a) 
dfi (b, R„ f ty/dn = BUb) i- K (b) e ‘'” S * (3-2b) 

(** , f, - cyr) + z C (r) e imB * 0-20 

where the coefficients are found, using Eq. (3-1) as input, from 

(3-3a) 
(3-3b) 


(3- 3c) 


A* m JL 

M "> ZTT 


**/» 


dp (*«, ^,0) ^ -inB0 d( £ 
<3 J? 


zr/s 

fl," - A. f J 4 . 

m l J~n ^ 9 


r n = JL 

" 27T 


,27r/a 


a^o('z d> T > 0) ^-c*rf30 


c/0 


£r~P( * 1r-Pn * ZfP? 

(3-lb) 


(3-lc) 
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The above forms apply to surfaces S A , Sb and S c respectively. For later 
reference we also quote here the analogous expansions for the pressure itself: 


where 


j>)= A a (r) + 2 Z A m fr) e. 

srr*/ 

Z (st) e c '” 3 ^ 

/n*/ 

■*■,$) = C.(r) + z&.£ C„C-h) e*'” 6 * ■ 

*?=/ 


d m (r)- Sr / ft) e ^ft 

o 

xrrfa 

3 M (r) = c/ft 

,2rr/a 

c »,(r) = jz- J f>(zj,+- y ft) e i ' n3< ^ c/ft 


(3-4a) 


(3-4b) 


(3-4c) 


(3-5a) 


(3-5b) 


(3-5c) 


Note that the subscript m on the coefficients is the Fourier index, while the 
absence or presence of the superscript n indicates a coefficient of either the 
pressure itself or its normal derivative. 


Equations (3-2) and (3-4) are written in a blade-fixed frame. The 
transformation to coordinates ( Z , ■F’ , <p ) which translate with the hub but do 
not rotate is simply 

2-2 -F - -f- ft * ft - Jit (3-6) 


for a propeller rotating in the -0 direction (as assumed by NASPROP-E) with 
constant angular velocity -/2 . When applied to Eq. (3-2) the result is 


&p(z« , ft, t)/#» 
d<p( z,# s ,ft, /)/ a* 

ft, £)/ 


a:<F) + Z AZ&)e im * r ** /l * ) (3 - 7a) 

Z a£ w 

c:m c;&e‘ m3 ^ /2 *{ 3.70 

^>f—/ 
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and similarly for Eq. (3-4). This clearly shows that the disturbance field on 
S, and by inference the entire field external to S, consists of a steady 
component plus a time variation at the Blade Passage Frequency , <3/2 , and its 
harmonics . 

Transformation from Mean Flow to Static Flow 


So far our attention has focussed on the boundary conditions to be 
applied to the outer flow at the interface S. We now consider the equations to 
be satisfied in this region, where nonlinear effects are assumed negligible. 

The linearized equations governing the outer flow including the effects of a 
subsonic axial Mach number, M, can be reduced to a single equation for the 
pressure , 

/ ) f>n * * - O (3-8) 

where /#* = l-M^, and subscripts indicate differentiation. This relation, 
written in the same hub-fixed coordinates as Eq. (3-7), is the convective form 
of the wave equation. As noted earlier, the BIE program is designed to work for 
the equations of classical acoustics with no mean flow. Setting M = 0 in Eq. 
(3-8), and further assuming a harmonic time dependence with frequency ^ = m 
in accordance with Eq. (3-7), we get 

-fti * f'C-Tf-)- + r -‘f> u y - O (3-9) 

which is the classic Helmholtz equation solved by BIE, k m being the wavenumber 
of the mth harmonic, m , 

Consider the following transformation of both independent and 
dependent variables (Refs. 10, 11): 

r * -F <p ** & 2-2 j/3 

exp (- i M k^zj/3) 

Applying this to Eq. (3-8) and cancelling common factors yields 

■fin- * + ( k ’»l/ 3 ) i fi - ° 


(3- 10a) 
(3-1 Ob ) 

(3-11) 
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Hence if we identify 



(3-10c) 


as part of the transformation, we recover the classic Helmholtz equation, (3-9). 
Thus we have transformed from a convective flow equation in (2 , -F , <f> » -fi* ) 
space to an analogous static flow equation in ( 2 , T , <p ,-p) space. An 
integral formulation can be used to solve the latter problem directly as 
follows . 


Briefly, the solution to Eq. (3-11) in the outer flow in Fig. 2 is 
uniquely determined by conditions on S and at infinity (outgoing waves only). 
Let denote the operator on the left-hand side, and G denote the Green's 
function satisfying the analogous inhomogeneous equation, 

oC G * >*J ~ (3-12) 

where ef is the Dirac delta function and x Q and x are generic symbols for the 
vector coordinates of the source and field points in ( 2. ,-?”,$>) space, 
respectively. As derived in any standard text on mathematical physics, e.g. 
Baker and Copson (Ref. 12), the solution to Eq. (3-11) can then be expressed 
in the following integral form, 



where the integral is over the whole surface S(x 0 ) and n 0 denotes the outward 
normal at the source point x Q on S. This form assumes only outgoing waves at 
infinity, and that there are no sources external to S. 


If a Green's function G can be found which vanishes on S, Eq. (3-13) 
reduces to a simple quadrature over p. This corresponds to replacing S by a 
dipole distribution. Conversely, if dG jdn 0 can be made to vanish on S, the 
solution reduces to a quadrature over 3f>/dn 0 , corresponding to a monopole 
representation on S. In general an analytical representation for G with either 
of these convenient properties is not possible for an arbitrary surface S, 
including the finite length cylinder used here. In this case it is 
convenient to use the free-space Green's function; either p or df>/dn 0 is 
assumed known on S, and Eq. (3-13) in the limit where x approaches the surface 
provides an integral equation for the remaining unknown. The inversion of this 
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integral equation is the function served by the BIE method described in Section 
4. The choice of whether to regard p or djo/dri,, as given in the BIE solution is 
here dictated by how easily the solved quantity could be incorporated in the 
outer boundary condition of NASPROP-E. From the discussion toward the end of 
Section 2, it is clear that df>/dn would be difficult to enforce as a boundary 
condition on the inner flow. So we have chosen to have djo/dn specified by the 
finite difference solution, as anticipated by Eqs. (3-1) through (3-7). 


It remains to be seen how the transformation in Eq. (3-10) affects 
the boundary conditions. First we note that the radius of the cylindrical 
interface S is unchanged by (3-10a), but its length is increased by the 
factor ' . The wavenumber of each harmonic of the acoustic field is also 
increased by this factor. Further, BIE will need as input not the dp /dr) 
specified in Eq. (3-7), but rather df>/dn in the transformed space. 
Differentiation of Eq. (3-10b) yields the following relation between the two: 


dp _ fdp d* . Mk„ dz 1 - 2 (3-14) 

JW ~ Ldn d» 1 fi, f J e r 


Thus we see that in order to specify dp/dh in the transformed space, in 
general both the pressure and its normal derivative must be supplied in 
physical space. Note that on S& and Sc, 


while on Sg, 


dn jSn = dzjdz - fi 


d/o j 

' a £ 

dr'jdr = / 

(3-14) specializes to 

■ZS- 

dn 

= 

/>!? - 

dp 

9PT 

= 


dp 

dn 

= 



dz[dn -?/ 
dzjdrf - 0 


; „ 7 TTST* ^ 


/> 


f]e ~1P 


iMk„ 


(3-lSa) 


(3-15b) 


(3-15c) 
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on surfaces S A , Sb and Sc respectively. Had an arbitrary curvilinear interface 
been used, the evaluation of Eq. ( 3— I 1 * ) would not have been so simple; it was 
for this reason that the cylindrical shape was chosen. It follows from Eq. 
(3-15) that the Fourier coefficients of df/an are related to those for p and 
dfildn ( C f. Eqs. (3-2) through (3-7)) by: 


i A1 


fop). I/3AZ(+-) + i^AJr)] e~ Z 3 ' 


(3-16a) 


LMk» 


dZ(z) = 3Z(z) e ‘ /3 s 


(3-16b) 


1/ik* 


czm =■ [/S eZM - i ££=■ c„ (+)] e ’ /3X 


(3-16c) 


which again apply to S A , Sb and Sc respectively. Note that when evaluating the 
above for m = 0, k m = 0. 


It is the Fourier coefficients defined by Eq. (3-16) which are 
actually passed to the BIE package as a representation of /&* over S. BIE 
then inverts the integral equation derived from (3-13) » as discussed in Section 
4, for the transformed pressure p on S. Actually, BIE must do an independent 
calculation for each harmonic, as the constant k m in Eq. (3-13) will be 
different for each value of m. It returns a set of Fourier coefficients A m , Bn, 
and C m representing the variation of p over S A , Sb and Sc in the transformed 
plane. These coefficients are converted to their counterparts in physical, 
blade-fixed coordinates through Eq. (3-10b), giving: 




. £ 


A„(rl - A m (?) 

e /* 


(3-17a) 


, l M kr>i 



B m (z) = S m (S) 

e /* 

£ 

(3-17b) 


. i Mkn, 



C„(-r) = S m Cf>) 

e A 

Z d 

(3-17c) 


When substituted into Eq. (3-4) these define a new pressure distribution over 
the surface S. 
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Update to Boundary Conditions 


The new pressure distribution found by BIE in general will not match 

that used as the boundary condition on the nonlinear flow equations solved by 

NASPROP-E. The disparity is used to update the boundary conditions used for the 

q 

next time step as follows. Let p£ denote the pressure distribution on S used as 

the boundary condition on the inner flow at the qth time step. This solution is 

used as input to the BIE package, as per the procedure described above, to 

determine the pressure over the same surface based on the outer linearized 
q 

equations, p Q say. The boundary conditions on NASPROP-E for the next time step 
are then updated as follows: 



ft? + <** (fto* 'ft*) 


= cLtf* + (/-eijflf 


(3-18) 


The parameter©^ f where 0<oC$/ , has been introduced to allow for 
underrelaxation of the iterations. This was found necessary for convergence of 
a similar iterative process used at Calspan in our studies of adaptive-wall 
wind tunnels (Ref. 13). As noted at the end of Section 2, p is not one of the 
primary variables used by NASPROP-E. Therefore, on S^ and S 5 the new values 
from Eq. (3-18) are used to recompute the local total energy at the boundary 
grid points using Eq. (2-4). On surface Sc, the Method of Characteristics 
(MOC) is used to update the boundary values for density, velocity, and energy 
using the new pressures. 


Eventually, the above iterative process should converge in the sense 

q q 

that the discrepancy between Pi and p 0 grows progressively smaller. Convergence 
is evaluated in CALPROP by testing whether 

[jr T Z (■??-■&)*] V e « ' (3-19) 

where the sum is over the total number, N? = JMAX+KMAX-1 , of grid points on S, 
and & is set by the user. When Eq. (3-19) is satisfied, the composite solution 
provided by NASPROP-E and BIE should consistently satisfy the blade boundary 
conditions, the free-field radiation conditions, and match p and across 

the interface S. 


The above procedure for interfacing between NASPROP-E and BIE is 
summarized briefly in Section 5 as a means of introducing the new routines in 
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rn 


the combined code CALPROP. This recapitulation makes clearer the sequence in 
which each of the operations is performed. 


After convergence the acoustic near-field can be evaluated by 
interpolation in the finite-difference grid. The acoustic far-field is easily 
evaluated from Eq. ( 3 - 13 ), which now becomes a simple quadrature over the known 
and dfi/dn distributions. This is also done in BIE, which is described in 
more detail in the next section. 
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Section 4 

OUTER LINEARIZED ACOUSTIC MODEL 


Review of the BIE Method 


A solution technique utilizing a discretized surface integral is 
commonly referred to as a Boundary Integral Equation (BIE) method. The 
outstanding feature of this method is that the dimension of the problem is 
reduced by one, since only the boundary (surface) of the body is dealt with. 
Thus, the BIE method is particularly suited for field problems involving an 
infinite domain or a finite domain in which the field variables are needed at 
only a few points. The BIE has been used successfully in solving problems 
governed by linear elliptic, parabolic, and hyperbolic partial differential 
equations (Ref. 14). Recently, Shaw (Ref. 15) has given an excellent review of 
the application of the BIE to wave problems. 

In principle, the BIE formulation can be derived from Green's theorem 
(Ref. 12, 16, 17) which relates the surface integral over S to the volume 
integral over V bounded by S for any two functions p and f/ sufficiently smooth 
and nonsingular in V: 

J(fj% ~ v %) (4 - 1) 

S If 

where n is the outward normal of S (i.e., n is directed out of the region V). 

We apply this identity to three-dimensional acoustic problems, e.g., 
sound radiation from a vibrating body of surface area S 0 , using the freespace 
Green's function ZXp (~L k £)/& = f* and the acoustic pressure p, where R is the 
distance between any two points P and Q in V, and k is the wave number. This 
problem is classified as an exterior problem governed by Helmholtz's equation 
where all physical quantities are defined outside of the body S 0 . The volume V 
is the space bounded by S Q and the surface K at infinity. The surface at 
infinity can be represented by a sphere of infinite radius. The acoustic 
pressure p is smooth and nonsingular, but f/ is singular when P = Q. To remove 
this singularity, we exclude a portion of V in a small sphere o' of radius £ 
surrounding P as shown in Fig. 3* Now, since P ^ Q, both f and p satisfy the 
Helmholtz equation V Z f + k 2 yj=0 and V p) + k z p> ~ O in V, the volume 
integral in Eq. (4-1) is 0: 
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Figure 3 RADIATION FROM A SURFACE S 0 INTO AN EXTERIOR REGION V 



(4-2) 



S 0 + <r+Z' 




o. 


Taking the limit as € tends to zero for the integral over o' 



) c/S 


— lim 

-O 


IJfrJ- 


ikP 


-ikfi 

e ap 


P an 



c/S = 4 Ttf> / P) 

(4-3) 


which results from using a spherical coordinate system centered at P and the 
fact that d/dn * - d/dP where p(P) is the acoustic pressure at P. The 
surface integral over the portion of S at infinity, i.e.,Z , is zero because 
of the Sommerfeld radiation condition. Defining a new normal V = -n directed 
outward from the body (i.e., into V), then Eq. (4-2) becomes 



= (P). 


(4-4) 


Now consider Q as a source point on the surface S Q and allow the field 
point P to approach S Q . Providing that the surface S Q has a unique tangent at 
Q, the singularity that occurs as Q = P may be removed by excluding the portion 
of V in a small hemisphere cf around P. Following the same procedure used in 
deriving Eq. (4-4), the following result is obtained: 

„ -iMQ -ikQ 7 

f[f b (nr) ~ V %]* - ^ rp> (4 - 5) 

Equation (4-5) is called the Surface Helmholtz Integral Equation. 


If the field point P is inside the body, then ip and p are always 
nonsingular in V. Therefore Eq. (4-1) yields, for this case, 
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which is called the Interior Helmholtz Integral Equation. Equations (4-4) 
through (4-6) can be written in a compact form: 


fb>k(g) - 
* 

where 

C(P)~H-T( for the field points outside the body, 

- ZfC for the field points on the surface of the body, 
= O for field points inside the body. 


As mentioned before, for field points on the surface S Q , C(P) - ZTt 
only when the surface has a unique tangent at P. A more general expression for 
C(P) is necessary if there is no unique tangent plane at P, i.e., when P is at a 
corner or edge of S Q . For this more general expression, note first that the 
counterpart of Eq. (4-7) for a function U satisfying Laplace's equation inside 
the region bounded by S 0 has the form 


f[w 3 S - u h(i)] dd ° = <**> d(p) 

4 


(4-8) 


in which (1/R) is the counterpart of exp(-ikR)/R for the Laplace operator, and 
C°(P) is the result of the limit of the integrals over a small "bubble -like" 
region O' with "base" centered at P similar to the hemisphere mentioned 
earlier. However for P at an edge or corner, the bubble-like region is not a 
hemisphere and the base is not flat in the limit. The required result can be 
obtained from Eq. (4-8) using U = 1 which gives 

c ° (p) --fk C4 - 9) 

which has the interpretation of the inner solid angle at P. Since the limit 
processes depend only on the order of singularity of \jj which is the same for 
both the Laplace and Helmholtz operators, and since C(P) in Eq. (4-7) may be 
obtained in the limit as P approaches S 0 from outside S Q , C(P) has the 
interpretation of the outer solid angle at P so that 
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C(P)+ C°CP) = VfC. 


(4-10) 


Thus, C(P) -tyfV - C°Cp) and Eq. (4-5) can be rewritten: 



Equation (4-11) is a general Surface Hemlholtz Integral Equation applicable to 
any arbitrary surface. 

Axisymmetric Formulation 


When the surface is a body of revolution (Fig. 4) then 

J>(Q)d<fi dL 0 (Q) 

where dL 0 is the differential length of the generator of the body. We now 
rename our normal ^ to n to conform with convention in acoustics. Further, if 
we expand p and djojdn on S 0 into Fourier series in the <p direction : 

/>7 



&-{Q) =LA" e imS * 

where B (the blade number) is a constant, then Eq. (4-11) may be written for a 

n 

single harmonic, A m and A m , as 



]fi dl. = CCP)A„CPJ 

(4-12) 


A //->»_ A 

where A^kP) — is the pressure at a field point P. The axisymmetric 

formulation above allows us to decouple the integration of Eq. (4-11) into the 
evaluation of two line integrals in the (L 0y <f>) coordinate directions, Eq. 
(4-12). Further, the integrals in the <fi direction are independent of the 
boundary data, A m and A m , and can be evaluated without the use of nodes and 
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elements. Thus, all element discretization can be confined to the generator of 
the axisymmetric body, i.e. the elements are now line segments. The boundary 
data need only be specified at the nodes on the generator. 


Numerical Implementation 

Because the exterior solution reduces to a discretization of the 
generator, only one -dimensional elements (line segments) are needed. The same 
quadratic interpolation is used for both the geometry and the function variation 
on each segment (isoparametric elements). The coordinates of any point 

on the generator (see Fig. 4) can be represented as: 

/>&-£**#)/* (4 - 13) 

z(€) - £ t Njt) ( 4 - 14 ) 

where * r ® the shape (interpolating) functions in the mapped domain 

(-/^dt^/ ) corresponding to node oc (each element is a line segment consisting 
of three nodes) and are the physical coordinates of node 06 . It is 

implicit here that all elements are mapped into a "clean" -1 to 1 region. All 
integrations are carried out in this new region. The Nd(&) are 

A 4 + ^4* , 

( 4 - 15 ) 

/v, ({,) -/-<?* , 

where £ is the coordinate in the -1 to 1 mapped space. 

n 

For an element k, A m and A m are represented using the same shape 
functions, Eq. (4-15): 

Amk (O - £ H. (i) ^ntkcL 

flC*/ 

( 4 - 16 ) 

AL (i) - £ (4) 

where A m UeL and A^hot. ar ® the pressure and normal derivative of the pressure, 
respectively, at node oC on element k for harmonic m. Substituting Eq. (4-16) 
into Eq. (4-12) yields: 
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(4-17) 


/ (Z) Hl (£)/> (t) 4 rt) 


-f 


-t a-^J'k^ rt)AU (£)/> (<*) s k {&) ctf]- cm/tjP ) , 


where 





ZfC 

a 

S>n 


_zrc 


W e , K„ (£)-_ /> « im * <*V> , 


(4-18) 


K is the total number of elements on the generator and is the Jacobian of 

the mapping from (jO } 2) to £ , given by: 


= [(j$) + (si) 

In the same manner, the constant C(P) in Eq. (4-17) is given by: 

c(f>) = * £ fx^(<t)f(iy k c^)dd , (4-i9) 




where 


*«(*)=/"§* (*)# 


(4-20) 


If we put P on the surface at a node which we identify as j 
(j = 1 to L, where L is the total number of elements) then A m (P) = A m j, and Eq. 
(4-17) may be written as 


K 

l, 


J 

L 

dL~f 


■* mUaL 


a. 


'Hi 


- A -j[* 


r., */ 


■ ~ 7 


J 


-// 

A=/ oL =f 


*n,/cct VmkJ > (4-21) 


«Z y - £ Al (4.22, 

-/ 

b:„j = £ f'K L (<t)K (t (4-23) 
w 

(4-24) 
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In the above expressions j denotes the global node number where P is located and 
k ,cc denote the element number and local node number, respectively, of the 
location of Q. If we let £ be the global node number for Q corresponding to a 
given k,o6 combination, Eq. (4-21) may be written as 

£ 4~! 4 ^ - J-jl* *£ ^kj7 - £ 4-4 > ( 4 - 25 ) 

or 

£ >0 - £ 4.* -O ( 4 - 26) 

where 

A' 

a ~f/ " ' A *£ 

where ^0/ is the Kronecker delta function. Equation (4-26) may be written in 

1/ 

matrix form as: 


[ Q n] { A ~} = [ b n,]( A Z.} ( 4 - 27 ) 

If the normal derivative of the acoustic pressure, ^A^j, is specified on the 
surface, then from Eq. (4-27) we may calculate the acoustic pressure {A m }. 

Calculation of Farfield Sound Pressure 


After A m has been determined we may determine the sound pressure at an 
exterior point P using Eq. (4-12), or its numerical counterpart, Eq. (4-17), 
with C(P) =l i-TC . This is a simple quadrature process since (now) both A m and A^ 
are known on the surface. 

The Interior Eigenfrequency Problem 

It has long been recognized (Ref. 18) that integral equations of the 
type given in Eq. (4-12) will not yield a unique solution at certain 
characteristic frequencies associated with the shape of S 0 . These 
characteristic frequencies are the eigenfrequencies of the interior of S 0 with 
an "altered" boundary condition. That is, if we are solving Eq. (4-12) by 
specifying dfjdn on S Q> then the characteristic frequencies where the nonuniqueness 
problem arises are the eigenfrequencies of the interior of S 0 determined by 
setting p = 0 on S Q . For a cylinder these frequencies are (Ref. 17, p. 429): 
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(4-28) 


***** - [(»*A)*+ (*»,«/& )*■]* 


where: n,<* = 1, 2, 3...; m = 0, 1, 2 ...; R and L are the radius and length of 

the cylinder, respectively; and Y mdL are the zeroes of the mth order Bessel 
function . 


The modal density (average number of characteristic frequencies per 
unit frequency) is proportional to k3. Even at moderate values of k the 
eigenfrequencies are closely spaced. The following table shows the 
eigenfrequencies between 19^ kR$ ^ 22 for the CALPROP interface (L/Rs = 1.21) 
for the m = 0 circumferential mode. 


oi 


^006 n^S 


2 

4 

6 

5 

3 

6 
1 
7 
2 
5 

4 
7 


7 

6 

3 
5 

7 

4 

8 


19.00 

19.54 

19.68 

19.79 

20.13 

20.84 

20.91 


1 

8 

6 

7 

2 


21.37 
21.49 
21.58 
21.66 
21 .84 


It has been observed (Refs. 5, 19 and 20) that symmetric boundary data 
will select only certain of the eigenfrequencies. Thus, only the eigen- 
frequencies above would result in nonuniqueness for symmetric (m = 0 ) ^0/^/7 
boundary data; i.e., the m = 1, 2 ... eigenfrequencies would not produce 
noriuniqueness . This is due to the interplay between the coefficient matrix and 
the right-hand-side vector (Ref. 19) in Eq. (4-27). In the present problem 
there is symmetry in the circumferential direction due to the blade spacing. 

Thus for eight blades only the ra = 8 modes will cause nonuniqueness problems. 
However, a calculation using Eq. (4-28) yielded 32 eigenfrequencies in the range 
99^kRg^101 for the CALPROP interface. 
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Attempts to Circumvent the Eigenfrequency Problem 


There is no general and systematic way to handle the nonuniqueness 
problem. This fact is evidenced by recent papers (e.g. Ref. 19) proposing new 
ways to achieve a unique solution to Eq. (4-12). However, the most common 
method (Refs. 18, 21) to handle nonuniqueness is called CHIEF. In CHIEF 
a least-squares solution is found by overdetermining the system, Eq. (4-27), 
with additional points interior to S 0 where p = 0. These additional equations 
are given by Eq. (4-6). 

A second method (Ref. 22) forms a new integral equation from a linear 
combination of Eq. (4-12) and its normal derivative. The normal derivative of 
Eq. (4-12) has a nonunique solution at certain characteristic frequencies, but 
these eigenfrequencies are not the same as those for Eq. (4-12). Thus, a linear 
combination of Eq. (4-12) and its derivative will have a unique solution for all 
values of frequency. This method has been used successfully recently with a 
combination finite element/boundary element solution approach (Ref. 10). 

Nonuniqueness of the solution to Eq. (4-12) is manifested as 
ill-conditioning in the algebraic system, Eq. (4-27). This ill-conditioning 
reduces the accuracy of the (unique) solution in a frequency band A k in the 
neighborhood of the eigenfrequencies. If it were not for this ill-conditioning 
over Ak , the nonuniqueness problem would be only of academic interest because 
it is always possible to avoid the eigenfrequencies by a proper choice of 
frequency. Thus, it is the ill-conditioning that results from nonuniqueness, . 
and not nonuniqueness per se , that produces inaccurate results in solving Eq. 
(4-27). 


The method used in Ref. 22 works when the eigenfrequencies of Eq. 
(4-12) and its normal derivative are sufficiently spaced so that the A k do not 
overlap. When this condition is met ill-conditioning will not occur. However, 
at higher frequency, e.g., kR > 20, it is expected that this method will not 
yield good results. It was this consideration that led us to pursue the CHIEF 
method over the method in Ref. 22. 
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Another method to overcome the fictitious eigenfrequency problem is 
the Exterior Overdetermination Method (EODM) (Ref. 19)* Based on a combination 
of the Surface and Exterior Helmholtz Integral Equations, EODM uses an iterative 
approach. The square system of equations obtained from the Surface Helmholtz 
Equation is overdetermined by additional equations formed by evaluation of the 
Exterior Helmholtz Integral Equation at selected points outside the body. An 
approximate impedance function is assumed on the surface of the body to initiate 
the iterations. The selection of the locations of the exterior points used in 
the iteration was found to be critical as a result of some numerical experiments 
we performed using EODM. If the points are too far removed from the body, two 
adjacent points yield almost the same equation, thus causing severe 
ill-conditioning. If the points are too close to the surface, it is not 
possible to have many points for the same reason. The optimal selection of the 
number of additional points and their locations corresponding to the desired 
rate of convergence remains unresolved. 


33 


Section 5 
PROGRAM CALPROP 

Sections 2 through 4 were devoted to the theory and algorithms used in 
the NASPROP-E and BIE codes and how they were intended to interact in the 
resulting combined program, CALPROP. This section is concerned more with the 
mechanics of how all this was implemented. Originally it was envisioned that 
CALPROP would consist of a new MAIN program and associated I/O which would 
essentially treat NASPROP-E and BIE as subroutine packages. However, it was 
soon realized that since NASPROP-E would comprise the majority of the source 
code, and already contained much of the required I/O structure, it was easiest 
to build the new program about its existing framework, with BIE incorporated as 
a package of subroutines. This approach has the added advantage of making 
CALPROP that much easier to run for previous users of NASPROP-E. Since 
NASPROP-E has been documented in Refs. 3 and 4, the emphasis here will be only 
on modifications to existing routines and descriptions of whatever new routines 
have been added. The source code images are stored on the NASA Lewis CRAY 
system as PDN = SRCLIB1 , and the executable object code as PDN = 0BJLIB1 ; both 
have ID = FSSPAN. The program takes up approximately 652K words of storage. 

The overall structure of program CALPROP is shown in Fig. 5. Although 
all subroutines in the progam are shown, this is not intended to be a detailed 
flow chart, but merely to display the major interconnections. The starting 
point was the NASPROP-E version running on NASA Lewis' CRAY computer, which was 
permitted to us for this work by Dr. L. Bober. Those routines outside the 
dashed line are native to NASPROP-E. Some of these required changes and are 
discussed first below. This is followed by a description of the new routines 
unique to CALPROP, i.e., those inside the dotted line in Fig. 5. Finally, a 
description of the BIE package is provided. 

Changes to NASPROP-E Routines 

The only existing NASPROP-E routines which required modification were 
AIR 3D, INPUT, IN IT I A and BNDRY as follows: 
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Figure 5. OVERALL STRUCTURE OF CALPROP 
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AIR 3D. This is actually the Main, or driving, program; the only 
change here is the addition of labelled COMMON /A CO UST/, which is the primary 
means of transferring data between the NASPROP-E and BIE routines. The 
variables, in the order which they appear, are defined below. 

Variable Name Description 


EPS 

Convergence criterion 

used in Eq. (3-19) 

RMS 

Left-hand side of Eq. 

(3-19) 

FSBETA 

(1_M2)V2 



Value of J grid index at upper left corner of 
cylindrical grid in Fig. 1 

JMAX+KMAX 

Highest value of harmonic index m, i.e., 0<m<MMAX, 
MMAXs* 5 

Integral number of time steps between boundary 
condition updates using BIE 

Number of observer positions at which the acoustic 
signal, SPL, is to be evaluated; NSIG ^20 

Relaxation factor on boundary conditions , oi , in Eq. 
(3-18) 

RSIG(IS) ,ZSIG(IS) Dimensionless r, z coordinates of the ISth observer 

position at which the acoustic signal, SPL, is to be 
evaluated; normalized by propeller diameter, D; 

1 $ IS S NSIG 


JA 


JMPKM 


MMAX 


NINT 


NSIG 


RLXBC 
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Variable Name 


Description 

SPL(MDUM,IS) The Sound Pressure Level of the mth harmonic at the ISth 

observer position, referenced to ; MDUM=m+1 , 

1 i MDUM £ MMAX+1 , 1 « IS « NSIG 


WAVK(MDUM) Dimensionless wavenumber of mth harmonic of acoustic 

field, k m , normalized by the reciprocal of propeller 
diameter 

DELPHI, PHASE Working arrays containing azimuthal grid spacing and 

phase information used in Fourier decomposition 

RTR(I) ,ZTR(I) Dimensionless r, z coordinates of Ith grid point on 

cylindrical interface, U 1^ JMAX+KMAX+1 


P(L,I),PN(L,I) Dimensionless pressure and normal pressure derivative, 

normalized by p* and p*, /D respectively; evaluated at the 
Ith grid point in the r", z plane and the Lth <p grid 
point. 


PC(MDUM,I), Dimensionless complex Fourier coefficients of the 

PNC(MDUM,I) pressure and its normal derivative, normalized by p^, and 

p^ /D respectively; correspond to the mth harmonic, 
m=MDUM-1, at the Ith grid point in the r, z plane. 


Labelled COMMON/ACOUST/ also was added to existing routines INPUT, INITIA, and 
BNDRY . 


INPUT. The principal change here was the addition of new input 
variables required by CALPROP . The original NASPROP-E would read in 4 input 
records from TAPE 5 for a start from scratch. In CALPROP, a fifth record is 
added to read in NINT, MMAX , NSIG and JA with a format 615, followed by a sixth 
record containing EPS and RLXBC with a format 6F10.5. Finally, the coordinate 
pairs ZSIG(IS), RSIG(IS) for IS = 1 ,NSIG are read in using as many input records 


37 



as needed with a 6F10.5 format, i.e., 3 pairs to a record. The values for NINT, 
MMAX , NSIG, JA, EPS and RLXBC are also written on the output file after the 
other input variables. 

INITIA. If NINT $ NMAX (cf. discussion of BNDRY below), this routine 
will now make initialization calls to new routines OUTER and BIE. 

BNDRY. There are two routines in NASPROP-E with similar names, BNDRY 
and BNDRY2. BNDRY2 updates the boundary conditions only on the blade surfaces, 
and was not changed in any way for CALPROP. BNDRY is responsible for updating 
the boundary conditions on the nacelle (no changes here), and on the outer 
surface of the difference grid, i.e. the interface surface Si As mentioned 
before, this surface was changed from the hemisphere/cylinder shape used in 
NASPROP-E to a finite length cylinder in coordinates. This 

modification only necessitated changes to the mesh generation program which 
produces the boundary-conforming grid. Hence to BNDRY, which works in the 
resulting (cf,)|,<^) coordinates, this particular change was transparent. 

As for updating the outer boundary conditions on S, BNDRY was altered 
so that how the boundary conditions are handled depends on the relative values 
of the three parameters NC, NMAX, and NINT. NC is the integral value of the 
current time step count, and is found in labelled COMMON/COUNT/. NMAX is an 
input which represents the maximum value of NC to be allowed on this run, and is 
found in labelled COMMON/BASE/; both are native to NASPROP-E. NINT is the 
interval, i.e., some integral number of time steps, at which BIE is called to 
update the boundary conditions. It is an input parameter communicated through 
labelled COMMON /ACOUST/ (see above), and is new to CALPROP. 

Three alternative logic paths are possible. First, if the user 
specifies his input values such that NINT > NMAX, BIE is never called, and BNDRY 
assumes a constant pressure of p^ over S for all time. Hence this represents a 
fallback option which will cause CALPROP to produce output identical to that 
from NASPROP-E. A second possibility occurs if NINT« NMAX, but NC is not an 
integral multiple of NINT. BIE will not get called on this step, and none of 
the boundary conditions on S are updated (though those on the nacelle will be). 
The third and final possibility is that NINT $ NMAX and NC is an integral 
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multiple of NINT as well. This triggers a call to BIE as well as several other 
new routines unique to CALPROP which are used to correct the boundary 
conditions. Thus, if NINT = 3, BIE will be used to correct the conditions only 
every third time step, while for NINT = 1 the conditions are updated on every 
step. The other new routines that are called in addition to BIE and the 
functions they serve are described below. 

The only other change made in BNDRY occurs when NC equals NMAX and is 
also a multiple of NINT. BNDRY then makes calls to NOISE (also described below) 
to compute the Sound Pressure Level, SPL, at those (RSIG, ZSIG) points inside 
the finite difference grid, and to the BIE package to get the SPL at points 
outside S. 


New Routines Unique to CALPROP 


These routines are called by BNDRY, directly or indirectly, when the 
step count NC is an integral multiple of NINT. They include OUTER, FRFIT, 
TRFORM, FREVAL, PR, and NOISE. We first give a brief description of each, 
followed by an explanation of how they interact. (The BIE package will be 
described separately.) 



OUTER 


FRFIT 


TRFORM 


Description 


Evaluates p and;, dfijdn in physical space at the outer 
cylindrical surface needed for the acoustic field; once 
the acoustic solution is obtained, this routine also 
updates the pressure on the cylindrical interface 


accordingly. 


Fits a Fourier series to p and 
of Fourier coefficients. 


dpjdn ; 


; output is the set 


Transforms Fourier coefficients of p and 




between 


the physical and transformed space, the latter 
corresponding to no mean flow. 
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Subprogram 


FREVAL 


PR 


NOISE 


Description 


Using the coefficients passed by BIE, this routine 
evaluates the Fourier series for p at the grid points 
on the cylindrical interface. 

A function subprogram which calculates static pressure 
at any grid point from the density, velocity 
components, and total energy stored in solution vector 
Q. 

Evaluates the Sound Pressure Level (SPL). for observer 
locations within the finite difference grid by 
interpolation . 


All of these subprograms, except PR, contain the labelled COMMON /ACOUST/, which 
is their primary means of communication. No attempt has yet been made to 
vectorize them, as their run times are much smaller than the already existing 
NASPROP-E routines. A few are called with a single input parameter, IGO, whose 
value determines which of several logic paths is to be followed on that call. 

How these routines interact is best explained by recapitulating in sequence the 
steps described in Section 3 for communicating between NASPROP-E and BIE; at 
each step the routine responsible for its execution will be identified. In this 
way we are able to provide a convenient summary of the procedure, and the order 
of presentation will parallel the logic flow in the code. 

A preliminary step is the initialization of various variables and 
arrays in COMMON /ACOUST/ and in BIE which can be calculated and stored once and 
for all. For this purpose NASPROP-E routine INITIA calls both routines OUTER 
and BIE with their input parameter IGO = 0. OUTER then calculates FSBETA, JMPKM 
and the wavenumber array WAVK; also arrays DELPHI and PHASE to be used 
subsequently in routines FRFIT and FREVAL. 


The most important function served by this call to OUTER is the 
filling of arrays RTR and ZTR. These are the transformed coordinates 'r, z of 
the outer grid points on S. They are obtained from NASPROP-E arrays Y and X, 
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corresponding to r and z, and Eq. (3- 10a). RTR(I) and ZTR(I) are 
one-dimensional arrays whose index I = 1 at the lower left corner of Fig. 1 (J = 

1 , K = KMAX ) , and increases as one moves clockwise around the circumference , 
attaining its maximum value at the lower right (J = JMAX, K = 1), where I = JMAX 
+ KMAX + 1. This is two more than the number of points actually on S, because 
at the corner points (J = JA, K = KMAX) and (J = JMAX, K = KMAX) the normal 
pressure gradient is double-valued. Hence the two points corresponding to I s 
JA and JA + 1 will have identical coordinates and pressures, but different 
values of djojdn • The same is true of the two corner points corresponding to 
I = JMAX+1 and JMAX+2. 

After exiting INITIA, and subsequently EIGEN (unchanged) , CALPROP 
enters the main loop in Fig. 5, which is over the time step counter, NC. Each 
time NC is an integral multiple of NINT, the following sequence is triggered 
from BNDRY: 

1 . OUTER is called with inp*ut parameter IGO = 1 . This signifies that arrays P 
and PN in labelled COMMON /ACOUST/ are to be filled with values for pressure and 
its normal derivative. Array P is filled by calling function PR for each outer 
grid point. PR simply solves Eq. (2-4) for p in terms of the stored values of 
density, velocity, and total energy. Array PN is filled by evaluating the 
appropriate form of Eq. (3-1) at each point, using second order accurate 
difference formulae for p^ , p^ and p^ . The values for I s JA and JA + 1 will 
differ because of the switch from Eq. (3-1a) to (3-1b), as will the values for 

I s JMAX+1 and JMAX+2 because of the switch from Eq. (3-1b) to (3-1e). 

2. FRFIT is called to do the Fourier decomposition. Equations (3-3) and 
(3-5) are approximated via the trapezoidal rule, using the discrete values of p 
and df> fdn stored in P and PN in Step 1 . PHASE holds the necessary array of 
complex phase factors in Eqs. (3-3) and (3-5), and DELPHI the uneven increments 
in the abscissa, . Both of these were computed in the initialization call to 
OUTER. The end result of FRFIT is the set of Fourier coefficients for pressure 
and its normal derivative, stored in arrays PC (MDUM,I) and PNC (MDUM,I) 
respectively. MDUM is an offset index for the mth harmonic, i.e., m = MDUM-1, 

1 ^MDUM^ MMAX+1 ; again, I refers to a particular point on S. 
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3. TRFORM is called with the input parameter IGO = 1. This specifies that 
the routine is to compute the Fourier coefficients of the normal pressure 
gradient in the transformed variables by applying Eq. (3-16) to PC and PNC in 
physical space. The new transformed values are stored in array PNC, overriding 
those computed in Step 2. 

4. BIE (discussed below) is called with input parameter IGO = 1. This 
signifies that the outer solution is to be obtained by inverting the appropriate 
integral equation, using the data in PNC as the inner boundary condition. BIE 
returns the solution for the pressure on S as a set of Fourier coefficients at 
the same grid points. The coefficients are stored in array PC, overriding those 
calculated in Step 2. Note again that BIE must essentially solve a separate 
problem for each of the harmonics, m=0,1..MMAX. 

5. TRFORM is called again, this time with IGO = 2. This indicates that the 
Fourier coefficients returned by BIE, which are in the transformed variables, 
are to be converted back to physical space. TRFORM does so by applying Eq. 
(3-17) to array PC. The new values override those returned by Step 4 as they 
are computed. 

6. FREVAL is called. The coefficients PC returned by Step 5 are those 
appropriate to Eq. (3-4) in physical space. FREVAL simply evaluates the 
Fourier series at all grid points on S. The resulting pressures are stored in 
array P, overriding those computed on Step 1. 

7. OUTER is called again, this time with IGO = 2. Its job now is to use the 
newly estimated pressures from Step 6 to update the outer boundary conditions to 
be used on the next pass through NASPROP-E. First the new pressure values are 
relaxed using Eq. (3-18). The resulting values are used in Eq. (2-4) to compute 
new values for the total energy on S a and Sg, assuming that the velocity 
components retain their free-stream values. The new energy values are stored in 
solution vector Q. 

8. Control is returned to BNDRY, where the relaxed pressure values on the 
downstream surface S c are used with the Method of Characteristics (Ref. 3) to 
calculate updated values for the density, velocity and energy. BNDRY also 
computes and prints out the r.m.s. change in pressure, as defined by Eq. (3-19). 
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9. If NC = NMAX, routine NOISE is called from BNDRY to compute the SPL at 
those observer locations, read in as coordinate pairs (RSIG, ZSIG), which lie 
within the finite difference grid. First, the quadrilateral cell in the 
computational (<£ , >^ ) plane containing the observation point is identified. 
Then bilinear interpolation is used to get the pressure, PP(L), at each of the 
corresponding azimuthal locations, PHI(L), L = 1...LMAX. A Fourier 
decomposition is then performed on this data, again using a trapezoidal 
approximation like that in FRFIT, to yield the complex amplitude for the mth 
harmonic, stored as CP. Finally, the SPL is computed from 


SPL = 20 to 9lo /CP/. 

Since the pressures, and henoe CP, have all been normalized by p^, , the above 
SPL is a relative scale referenced to p^ . To put it on an absolute basis once 
p m and a reference pressure, p re f, have been specified, one has merely to add 
the conatant value, 20 logio (P,*, /Pref^* NOISE concludes by writing the 
coordinates of the observation point and the SPL predicted there for each 
harmonic. For observer positions outside the finite difference grid, the 
acoustic signal is computed with a call from BNDRY to BIE with input parameter 
IGO = 2. 

Modified Mesh for the BIE 


The BIE method involves integrations of kernels (e.g., exp(-ikR)/R) 
which are of oscillatory nature. The oscillations increase with frequency. An 
estimate of the length for an element to perform satisfactorily accurate 
integrations can be determined. It is assumed that the kernel exp(-ikR)/R can 
be integrated accurately, along with the shape functions, in an element of 
length equal to one-half the wavelength using five Gauss points per element. 
Thus, an element length of Tt/k is adequate. This is considered to be a 
conservative estimate of the element length. 

The high frequencies (kRg £ 100) involved require much smaller 
element size than some of the elements in the NASPROP-E mesh. Since the BIE 
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code requires quadratic elements and the NASPROP-E mesh did not conform to such 
a mesh, it was decided to generate a new mesh based on the NASPROP-E mesh. 

The NASPROP-E mesh is optimized for variations in boundary data, the 
nodes being crowded in regions of high gradients. Hence, it was decided to 
modify this optimal mesh rather than to generate a new mesh without 
consideration of boundary data. 

It was also observed that the NASPROP-E mesh contains two distinct 
groups of two-noded elements, the "wide" and "narrow" elements. The elements 
longer than SMALEL=0.5 (SMALEL may be changed as noted below) are considered to 
be wide elements and the rest as narrow elements. The NASPROP-E mesh also 
incorporates "dual nodes" at the corners of the cylinder. The BIE requires 
adjacent elements not to differ very much in length. Considering all these 
factors the following mesh generation scheme was adopted (see Fig. 6). 

One of the dual nodes was dropped at each corner. One node was 
introduced in the narrow elements and an odd number of nodes (NPTS) was 
introduced in the wide elements. NPTS is required to be an odd number so that 
an integral number of quadratic elements is generated in any wide element. 

One node was introduced in the narrow elements to make them quadratic elements. 
In this way, all the elements involved have their mid-nodes in the middle of 
the element. The narrow elements could not be combined (without addition of 
more nodes) into quadratic elements since an element might have been generated 
at the point of transition from the narrow to the wide elements, in which case 
the mid-node might have been very much off-center. This would cause singularity 
problems in the shape functions. 

Extra elements of the same size as the last narrow element were 
introduced in the nacelle to close the grid , and dj&jdn was prescribed to be zero 
on these elements. The given boundary data were approximated by a piecewise 
spline fit and the boundary data for the newly introduced nodes were determined. 

Description of Routines in BIE 

The block labelled BIE in Fig. 5 actually represents a package of 
several subprograms used to solve the outer acoustic field. A schematic of how 
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they are interrelated is given in Fig. 7, and below is a brief description of 
the function served by each. 

Subprogram Description 

BIE Driver program for the rest of the package'. IGO = 0 

represents an initialization call; for IGO = 1 a new 
outer solution is obtained; for IGO = 2 the SPL at 
points outside the interface is computed based on the 
last solution. 

SRI Called only by BIE. Puts exterior coordinates 

(initially in arrays RTR & ZTR) into array XI. 

SR2 Called only by BIE. This is the new mesh generation 

routine. Generates new nodes between existent nodes. 
Drops dual nodes at corners. Establishes node 
correspondence via arrays C0RR1 and C0RR2 between the 
given mesh and the new generated mesh. Introduces extra 
nodes in nacelle. 

SR5 Called only by BIE. Calls SPLFIT and function FS. 

Interpolates the given boundary values of the given mesh 
to determine the same for all nodes of the new mesh. 

Also establishes the local-global node correspondence 
for the new mesh. 

HEAX Called only by BIE. Calls GAUSS, SHFUN , COEF, EXTR , 

SOLVER and BIVL. Solves for the unknown surface and 
field pressures using the BIE method. 

SR 6 Called only by BIE. Calculates the Sound Pressure Level. 

FILLPC Called only by BIE. Fills array PC, i.e., finds the 

Fourier pressure coefficients on the surface at the 
nodes of the given mesh. 
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Subprogram 


Description 


GAUSS 


Calculates abscissas and weight factors for numerical 
(Gaussian) quadrature. 


SHFUN 


Calculates and stores values of shape functions and 
their derivatives for quadrature. 


COEF 


EXTR 


Calculates and stores the matrix of coefficients and the 
known {i.Q.,dpjdn ) vector. 


Modifies the matrix of coefficients and the known vector 
for exterior problems. 


SOLVER 

BIVL 


SPLFIT 

FS 


Solves the system of linear algebraic equations. 

Prints nodal values of p and dpjdn and prints the values 
of pressure at exterior points. 


Interpolates 


a^a 


to get cubic spline coefficients. 


Uses spline coefficients from SPLFIT to determine dpjd?) 
at additional nodes. 


Description of Important Parameters in BIE 

The following are major parameters governing the operation of the BIE 

solution : 


Variable Description 

NG Number of Gauss points per element for integrations 

along the generator when both source and field points 
are on the surface. Maximum value is 10. 
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Variable 


Description 


NGG Number of Gauss points per element for integrations in 

the direction of the angle of revolution. Maximum value 
is 10. 

NPTS Number of nodes introduced in every wide element in 

NASPROP-E grid. Must be an odd number. 

NSS Number of subdivisions (elements) in TC radians in the 

direction of revolution. 


SMALEL A length such that elements of NASPROP-E grid longer 

than this length are classified as "wide" elements. A 
value of 0.5 is suggested. 


The first four parameters are related - all affect the numerical 
accuracy. We recommend setting NG and NGG to a value between 5 and 10. When 
this is done we have found that a conservative element (subdivision) length is 
given by n/k where k is the wavenumber. This criterion can be used to select 
NPTS and NSS such that high numerical accuracy is achieved. Using the ftjk 
criterion we determine NSS by: 


NSS 


one-half of circumference of cylinder 
It/k 


+ 1 , 


NSS = 5k + 1 , 


since Rs=5 is the radius of the cylinder. The 1 in the above equation insures 
that there is a single subdivision in ft radians in the direction of revolution 
when k=0 (i.e. the m=0 term of the Fourier series). We use the same ft/k 
criterion to determine the maximum element length along the generator. We do 
this by selecting the wide element with the maximum length ( » 1 ) , divide this 
length by ft/k and pick NPTS accordingly (noting that NPTS must be an odd 
number). The relationship between k, NPTS and the total number of nodes on the 
generator is given below (for SMALEL = 0.5). 
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k< 


NPTS 


TOTAL NUMBER OF NODES 


It 

1 

139 

2ft 

3 

161 

3 it 

5 

183 

4ft 

7 

205 

5ft 

9 

227 

6ft 

11 

249 

7ft 

13 

271 
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Section 6 
RESULTS 


We have not been able to get the new CALPROP code to run a case 
through to completion, and so have little in the way of numerical results to 
report. Basically, the code takes too long to execute. In the discussion 
below, we will first describe the sample case we have been working on, followed 
by a description of the difficulties we encountered. 

Figure 8 displays the input data as printed by the code for the case 
we worked with, corresponding to an eight-bladed SR-3 propeller at a nominal 
flight Mach number of 0.8. The angle of attack is 0° (neither NASPROP-E nor the 
BIE package can handle the asymmetric inflow case), and the advance ratio is 
3.06. The finite difference grid used was that generated for us by H. Huynh 
(cf. Section 2) based on the SR-3 blade geometry and the new cylindrical outer 
boundary. Referring to Fig. 1, the radius of the outer boundary was set at 
Rs =5 propeller diameters. The upstream and downstream faces of the cylinder 
were put at z u s -5 and zg s 1.05. In the coordinate system used by the outer 
solution, these transform to z u = -8.33 and z j = 1.75. The grid has JMAX = 45 
points in the £ direction, KMAX = 21 in the direction, and LMAX = 11 in the £ 
direction (not shown). The line of constant <£ intersecting the upper left-hand 
corner of the cylindrical interface is at J s JA = 7. The grid point 
coordinates of this mesh can be accessed through permanent dataset name OURMESH 
with ID = FSSPAN. Because this was to be a demonstration run, the maximum 
number of iterations (time steps) asked for was only NMAX =5. A Courant 
number of 20 was specified for calculating the step size used by the ADI 
difference scheme. 

The new input parameters unique to CALPROP are listed near the end in 
Fig. 8. NINT = 1 specifies that the BIE package is to be called on each time 
step to update the outer boundary conditions applied to the inner flow. 

MMAX s 2 specifies that the m = 0, 1 and 2 Fourier components are to be included 
in the acoustic calculation. These correspond to the steady-state component, 
the fundamental blade passage frequency, and its first harmonic, respectively, 
cf. Eq. (3-7). This is considered minimal; the program is dimensioned to 
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MACH NO.* 0.80 

AN6LE OF ATTACK* 0.00 DE6REES 

RATIO OF SPECIFIC HEATS=1.4 

PINF=0. 1000E+01 

RH0INF=0. 1000E+01 

ADVANCE RATIO =-0.30600000E+0l 

BETA 3/4* -58.70 DE6REES 

BLADE DIAMETER* 1.00 FEET 

ONEGA* -0.31 REVOLUTIONS PER SECOND 

NO. BLADES* 8 

PERIOD = 45.00000RADIANS 

JMAX=45 

KMAX=21 

LNAX*ii 

ITERATIONS* 5 

NETH0D*2 ( 1*EXPLICIT, 2*IHPLICIT) 

COURANT NO. *20. 00 

4TH ORDER SMOOTH C0NST=0.1000 

2ND ORDER DISSIP CQNST=0.2000 

DIFFERENCING NETH0D*4 !2*SEC0ND ORDER, 4*F0URTH ORDER! 

IREAD*0 

IHRIT-1 

IPTCH*1 

JLE=20 

JTE*35 

KTIP=12 

IBC* 0 

INETOT* 0 

NCHNGE* 9000 

ISRID* 1 

HINT = 1 

UMAX * 2 

NSIG * r 

JA * 7 

EPS * 1.000E-02 
RLXBC = 1.000 


Figure 8 INPUT DATA FOR DEMONSTRATION CASE 
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allow MMAX to be as large as 5. Values of the other input parameters are not 
critical to the following discussion. 

After first mating the NASPROP-E and BIE codes together with the new 
routines required for the matching (Section 5), it was found that the combined 
CALPROP code would not execute. What ensued was a period of several months 
during which a series of trial-and-error runs were made to "get the bugs out." 
These included an initial incongruity between NASPROP-E and the new fully 
cylindrical grid, subtle logic errors which passed muster at compile time but 
raised havoc during execution, and incompatiblities between NASPROP-E and BIE 
which became evident only when passing information between the two. Without 
going into detail, suffice it to say that this initial debugging period took 
many times longer than anticipated. 

In our defense, we should point out that NASPROP-E and BIE were 
developed independently for somewhat different purposes, and on different 
systems at that. Our efforts were also hampered by the geographical separation 
between ourselves and Prof. Seybert, the time needed to familiarize ourselves 
with the idiosyncrasies of the NASA/Lewis CRAY operating system, and the 
limitations imposed on us by having only remote access to the system. The 
limited number of telephone access ports relative to the apparently large 
community of users pretty much restricted our useful window of access to the 
CRAY to a few very early-morning hours. 

Our preliminary development efforts have now brought us to the point 
where the code runs, but a more serious obstacle has been encountered. This was 
discovered when a run bombed for having exceeded our self-imposed job time 
limit of 100. secs., without completing the first time step. As a point of 
reference, the original NASPROP-E code completed 10 time steps in approximately 
12 sec. using the same finite difference grid. We activated the TRACEBACK 
option on subsequent runs , and added internal diagnostic messages to let us know 
which of the new routines unique to CALPROP had been successfully exited. This 
pinpointed the problem to routine COEF, part of the BIE package supplied by 
Prof. Seybert, which was entered but never exited. (This was the first time the 
code had gotten far enough along to call COEF.) As a further check we also made 
a run in which the entire BIE subroutine package was removed and replaced with 
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a ’’dummy" BIE routine that simply filled the PC array with the complex constant 
(1.0, 1.0). Though the results were meaningless, the modified code made it 
through two complete time steps in a few seconds. This demonstrated that the 
other new routines in CALPROP (dashed outline in Fig. 5) ran in a reasonable 
amount of time. 

Our initial suspicion was that perhaps a DO LOOP parameter in COEF had 
been left undefined, and the code was looping ad infinitum. A careful 
inspection of the COEF structure turned up the following set of five nested DO 
LOOPS : 


DO 500 K = 1 ,M 
DO 100 IG * 1, MO 


DO 90 IR = 1 ,N 


DO 515 IJ = 1 ,NS 
DO 51 IIG = 1 1 ,NGG 


51 CONTINUE 
515 CONTINUE 


90 CONTINUE 
100 CONTINUE 


500 CONTINUE 
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A system timing routine was added to COEF on an ad hoc basis, which confirmed 
that the program was in this section of code when it ran out of time. The 
symbolic dump triggered by the TRACEBACK disclosed that M = 89, NG = 10, N = 

179, NS = 80 and NGG = 10. This means that the block of statements within 
the innermost loop would have to be executed approximately 1.3 x 10® times. The 
block contains on the order of 10^ floating-point operations, which is probably 
on the low side. Hence a total of 1.3 x 10^ operations are needed; assuming a 
sustained rate for the CRAY of 2 x 10? floating-point operations per second, 
this section of code alone would require in excess of 10 min. CPU time. 

The above estimate is just for one Fourier component at one time step. 
For the present demonstration case which includes only three Fourier components , 
and assuming that oh the order of 100 time steps would be needed to converge 
(typical of the original NASPROP-E code, cf. Ref. 3, p. 100), something in 
excess of 53 hrs. of CRAY time would be required to get one solution. 

Admittedly this is only a ballpark figure, but it is clearly unacceptable even 
for a research oriented code. 

COEF is not a peripheral routine which can be temporarily bypassed. 

It is at the heart of the BIE package, in that it calculates the coefficient 
matrix in the set of linear equations to which the original integral equation is 
transformed by the discretization. The elements of the coefficient matrix 
involve quadratures in two dimensions over the cylindrical interface. It is 
these quadratures, for each of the M elements on the generator of the cylinder, 
which are being done in this nested loop construct. Parameters such as NGG, NS, 
etc., refer to the number of Gauss quadrature points and subdivisions to be used 
in the integrations (cf. Section 5). 

We consulted Prof. Seybert regarding our problem, who felt in 
retrospect that perhaps he had been overly conservative in choosing the number 
of discretization elements and Gaussian quadrature points. As a result, the BIE 
package was restructured so that a different number of interface elements and 
azimuthal subdivisions was used for each frequency (the lower the wavenumber, 
the coarser the grid), and NG and NGG were both reduced from 10 to 5. 
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Another attempt was made with these changes to CALPROP, and this time 
the first call to COEF was successfully completed for the m = 0 (i.e., 
steady-state) Fourier component in only 15 sec. However, for the m = 1 
component (the fundamental blade passage frequency), the code again bogged down 
in COEF and ran out of time, taking 175 sec. just to complete the K = 6 pass 
through the outermost loop, out of the 129 passes required. On this call, 

M = 129, NG = 5, N = 259 , NS = 10, and NGG = 5, for a total of approximately 4 x 
107 passes through the innermost loop. This represents only a marginal 
improvement from the previous run, and does not augur well for what would be 
required if it ever got to the m = 2 harmonic. 

Thus in its present form the CALPROP code can be made to run in some 
sense, but it is economically unfeasible to use. This is not because of any 
known errors per se, nor is the problem amenable to a "quick fix" resolution. 

The way the code is structured, it simply generates such a volume of arithmetic 
for the outer solution that it overtaxes the capabilities of even a CRAY. 
Unfortunately, at the point this was discovered there was neither the time nor 
resources remaining to do any major restructuring of the BIE package. More will 
be said about this in the next section. 
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Section 7 

CONCLUSIONS AND RECOMMENDATIONS 


Since we were not able to get any converged solutions with CALPROP, no 
meaningful conclusions can be drawn regarding the aerodynamics and acoustics of 
advanced propfans per se. Nevertheless, we still feel strongly that such a 
hybrid numerical scheme is a viable approach for such predictions, and would 
like to offer some comments that may prove useful to future investigators. . 

This is the type of problem that would ordinarily form the basis for 
a long-term program of developmental research. Instead, because of time and 
manpower limitations, we were forced to turn to "off-the-shelf" software 
wherever possible. In retrospect we may have been overly optimistic in thinking 
we could easily marry two such disparate codes as NASPROP-E and BIE. Each had 
originally been developed for independent purposes, and on different systems. 
NASPROP-E was attractive because it was under no proprietary restraints, and the 
mechanics of how to fold in the complex SR-3 blade geometry had already been 
done for us. Similarly, BIE was chosen because it was readily available, and 
had already been used with a finite cylindrical boundary. 

But neither code had been designed to communicate with another. 

Hence there were many inconsistencies between the two, some documented and 3ome 
not, that had to be ironed out on an ad hoc basis. This and the logistical 
problems introduced by the long-distance collaboration with Prof. Seybert, and 
having only remote access to the CRAY system, diverted too much effort away from 
the real technical issues. Ideally, a program of this complexity should command 
sufficient resources that the inner and outer flow solution packages would be 
developed in concert, on the same system, with the task of communication between 
the two always kept uppermost in mind. 

Having said that, there remain two obstacles to further development of 
the CALPROP code, viz., its very large run times, and the potential 
ill-conditioning of the outer solution near the interior eigenfrequencies 
(Section 4). The root cause of the run time problem is routine COEF, which 
calculates the coefficient matrix in Eq. (4-27) anew at each time step. Yet 
there is no reason to do so, since the matrix depends only on the acoustic 
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frequency and the node placement, which remain fixed. It is only the right-hand 
side vector in Eq. (4-27), which contains the boundary conditions, that must be 
updated on each iteration. Hence it is recommended that BIE be restructured so 
that the matrix is calculated and stored once and for all at the beginning of a 
run, preferably on the initialization call from INITIA with IGO = 0. Since a 
separate matrix must be stored for each Fourier component, this will greatly 
increase storage requirements . But with the memory available on the CRAY it 
should prove feasible, and the savings in run time would be tremendous. 

Even with only a single calculation of the coefficient matrices, the 
timing statistics quoted in Section 6 suggest that the run time might still be 
large. Further reductions could be gained by reducing the number of nodes used 
on the interface, and hence the size of the matrices, in the outer solution. 

The BIE package presently uses the grid intersection nodes from NASPROP-E as a 
starting point, and then inserts extra nodes between them as needed, cf. Section 
5 and Fig. 6. This is done to avoid abrupt transitions from closely to widely 
spaced nodes, which causes problems in the BIE method. But even the modified 
node structure that results is hardly optimal from the BIE viewpoint. It is 
suggested that a better approach is to base the node placement for the outer 
flow primarily on the needs of that solution, more or less independently of the 
grid points used in the inner flow. Cubic splines could be easily used to 
transfer data from one set of points to the other as needed. Such a scheme 
should allow fewer nodes for the BIE solution, and reduced run times. 

These last two points, i.e., when should the BIE matrices be 
calculated and for what node distribution, are good examples of the problems 
that arose because BIE was originally developed with other applications in mind. 
It had previously been applied only to comparatively coarse grids, where the 
boundary conditions had a very simple variation that was prescribed a priori; 
hence no iterations were necessary. Run times were minimal, so that 
considerations of computational efficiency such as those above never arose. 

Assuming the time needed for the outer solution can be reduced 
sufficiently, that still leaves the question of nonuniqueness. That is, the 
matrix equation for the exterior Neumann (Dirichlet) problem is known to be 
strongly ill-conditioned near the eigenfrequencies of the complementary interior 
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Dirichlet (Neumann) problem. As discussed in Section 4 for the frequency range 

r*t 

of interest here, kRg ■£> 100, the eigenfrequencies are much too closely spaced 

to get solutions by simply "tweaking" k until it lies betwen adjacent 
eigenfrequencies. For the present program the CHIEF method (Refs. 18, 21) was 
chosen as the most promising scheme to overcome this difficulty. Test cases 
have been run for simple point sources within the same cylindrical interface 
used here, but to date accurate solutions have only been obtained up to 
kRg— 20. This is well below the range needed for even the fundamental blade 
passage frequency, which for the case cited in Section 6 is kiRg = 109.5. For 
this reason the additional routines required to implement CHIEF are not 
presently included in CALPROP. 


Professor Seybert's statements at the end of Section 4 
notwithstanding, the author believes that the method used by Meyer, et al. in 
Ref. 22 is a preferable means of removing the ill-conditioning. They solve a 
modified integral equation, formed as a linear combination of the original 
equation and its derivative normal to the interface. The new equation has been 
shown to possess a unique solution at all wavenumbers. Numerical results 
presented in Ref. 22 confirm this, even at conditions coincident with the 
interior eigenfrequencies. This method was successfully implemented by 
Baumeister and Horowitz (Ref. 10) in a similar hybrid scheme applied to turbofan 
inlet acoustics. A finite element solution was applied to the inner nonlinear 
flow. The boundary integral method of Ref. 22 was used for the linearized 
acoustic field. Assuming (conservatively) that the radius of their interface 
was at the outer surface of the nacelle, for their blade passage frequency 
kiRs-35. The computed results were obtained with only 59 boundary segments and 
show no evidence of ill-conditioning. Reasonable agreement with experimental 
data was also demonstrated. 

The nonuniqueness of the solution to Eq. (4-12) near the 
eigenfrequencies is not a manifestation of any physical phenomenon. It is a 
purely mathematical artifact introduced by the transformation of the governing 
equation from differential to integral form. The fact that the problem gets 
progressively worse with increasing frequency should perhaps suggest to us that 
we may be attacking it from the wrong direction. That is, for the regime of 
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interest here we may not need the full-blown BIE artillery, if only we 
are clever enough to somehow take advantage of the short wavelengths. 

Such a high frequency approximation has been used successfully in 
studying radiation from solid bodies, e.g., Ref. 23 and 24. In the high 
frequency limit the waves behave asymptotically like plane waves. This 
together with the condition of no flow through the wall leads to a simple 
algebraic relation between the acoustic pressure and the velocity normal to the 
surface, which in effect can serve as a boundary condition replacing the 
integral equation. 

Such an approximation was briefly considered here, until it was 
realized that this simple relation would not apply in our case because the 
surface is not solid. For a solid radiator in this limit, the normal to the 
surface is also normal to the acoustic wave fronts; but in the present problem 
there is no reason for this to be so, as waves can pass through the interface at 
any angle. This added degree of freedom destroys the simple algebraic relation. 
However, it may be that with enough thought a more general relation could be 
worked out. The possibility of such an approach at high frequencies is very 
attractive, as it not only removes any ill-conditioning, but obviates 
altogether the need to invert an integral equation. 

In summary, the lack of success in the present investigation should 
not be taken as evidence that the basic approach is unworkable. A hybrid scheme 
matching an inner nonlinear flow to an outer linearized field still appears 
ideally suited, to predicting the aerodynamic and acoustic fields of advanced 
turboprops. It is hoped that the comments offered here will help to further 
progress toward that goal. 
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NOMENCLATURE 


a speed of sound 

A m , C m dimensionless Fourier coefficients of pressure distribution over 
S At Sg, Sc; defined by Eq. (3-5) 

A m , §m> Cm dimensionless Fourier coefficients of pressure distribution in 
transformed variables; related to the above by Eq. (3-17) 

^m* ®m» Cm dimensionless Fourier coefficients of the normal pressure gradient 
over S A , Sgi Sc; defined by Eq. (3-3) 

** n *** n rJ n 

A m> Bjj, C m dimensionless Fourier coefficients of the normal pressure gradient 
in transformed variables; related to the above by Eq. (3-16) 

B number of blades 

D propeller diameter 

e total energy per unit volume normalized by p^ ; defined in Eq. 

(2-4) 

E, F, G flux vectors in the £ , directions, respectively; defined in 

Eq. (2-6) 

G Green's function defined in Eq. (3-12) 

H vector of undifferentiated source terms defined by Eq. (2-8) 

J Jacobian of the generalized coordinate transformation defined by 

Eq. (2-5) 

J, K, L integer grid indices in the , y[ , £ directions 

JMAX , KMAX , 

LMAX maximum values of the above - 
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M 

m 

n 

P 

a/ 

p 

Q 


dimensionless wavenumber of mth harmonic, normalized by D 
free -stream Mach number 
Fourier series (harmonic) index 

local outward normal to the interface surface S in Fig. 2 
dimensionless static pressure normalized by 

dimensionless static pressure in transformed variables; related to 
p by Eq. (3-10b) 

vector of dependent variables in finite difference solution, Eq. 
(2-3) 


*s 


dimensionless radius of surface Sg, normalized by D 


S A» S B» S C upstream, sidewall, and downstream faces, respectively, of 
cylindrical interface S, Fig. 2 


t 


dimensionless time normalized by D 



u, v, w 


dimensionless velocity components in the 2, ■f' t directions 
normalized by /VJ" 1 


U, V, W 


dimensionless contravariant velocity components in the £ , >7 , £ 
directions normalized by CL g 0 j^Y' » defined by Eq. (2-7) 


dimensionless cylindrical blade-fixed coordinates normalized by D 


dimensionless cylindrical coordinates translating with the 
propeller, but not rotating; related to the above by Eq. (3-6) 

Z } -r ,<p dimensionless cylindrical coordinates used in transformation to no 

mean flow; related to the above by Eq. (3- 10a) 


64 



z u > z d 


oc 


/ 

/ 

€ 

/> 

n 


n 

q 

T 


m 


i 


o 


o 


CO 


dimensionless axial coordinates of S A and Sc, Fig. 2 

relaxation factor used in updating the boundary conditions on S, 
Eq. (3-18) 

a/ 1-M 2 ' 

specific heat ratio 

Dirac delta function 

convergence criterion, Eq. (3-19) 

generalized boundary-conforming coordinates, Eq. (2-1) 
dimensionless density normalized by 

dimensionless angular velocity of the propeller, normalized by 

a./foVF) . 

Superscripts 

used to indicate Fourier coefficient of the normal pressure 
gradient, e.g. Aq 

time index in Eq. (3-18) 

indicates the transpose of a vector, e.g. Eq. (2-3) 

Subscripts 

pertaining to the mth harmonic of Blade-Passage Frequency 

pertaining to the inner flow solution 

pertaining to the outer flow solution 

denotes a source point on S in Eq. (3-13) 

evaluated at upstream infinity 
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